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ABSTRACT 

A stochastic model of turbulent mixing has been developed 
for a reactor in which mixing is represented by n-body fluid 
particle interactions (n = 2,3,.. .6). The model has been used to 
justify the assumption (made in previous investigations of the role 
of turbulent mixing on burner generated thermal nitric oxide and 
carbon monoxide emissions) that for a simple plug flow reactor, 
composition nonuniformities can be described by a Gaussian distri- 
bution function in the local fuel:air equivalence raito. Recent 
extensions of this stochastic model to include the combined effects 
of turbulent mixing and secondary air entrainment on thermal gen- 
eration of nitric oxide in gas turbine combustors are discussed. 
Finally, rate limited upper and lower bounds of the nitric oxide 
produced by thermal fixation of molecular nitrogen and oxidation of 
organically bound fuel nitrogen are estimated on the basis of the 
stochastic model for a plug flow burner; these are compared with 
experimental measurements obtained using a laboratory burner operated 
over a wide range of test conditions; good agreement is obtained. 





I. INTRODUCTION 

In most practical continuous flow burners, where the fuel and 
the air enter separately, primary combustion is generally observed 
to occur within a relatively confined, turbulent, mixing-controlled 
flame region. The product gases which issue from this primary zone 
are exceedingly nonuniform both in terms of their species composition 
and temperature and, in the absence of heat transfer to the surroundings, 
will only approach the composition of the corresponding adiabatic 
premixed system after turbulence and diffusive mixing have made the 
flow more uniform. 

There have been several studies carried out recently, notably 
by Fletcher, Heywood, and their co-workers, e.g. [1-3], which have 
attempted to examine the role of turbulent mixing within the product 
gases on burner generated nitric oxide and carbon monoxide formation 
using the plug-flow reactor concept proposed by Be^r and Lee [4]. 

The important elements of the method developed by these investigators 
for quantitative modelling of this particular aspect of burner 
emissions problems are most simply illustrated by the study reported 
by Pompei and Heywood [3]. These authors carried out an experimental 
investigation using a relatively simple cylindrical burner in which 
an atomized spray of kerosene fuel was injected through a nozzle, 
located on the burner axis, around which combustion air at atmospheric 
pressure was admitted through swirl vanes. By sampling the mean 
product gas composition over the cross-section of the burner at 
several downstream axial positions, it was established that the 
flow was essentially one-dimensional in the mean over the major 
portion of the burner length. Therefore, Pompei and Heywood modelled 
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the flow throughout the burner as a plug-flow by assuming that the 
fuel was rapidly dispersed throughout the burner cross-section and 
burned within the primary zone to produce product gases composed of 
a large number of fluid elementSj or eddies, each characterized by 
a local value of the fuel:air equivalence ratio, (Ji, and each having 
a scale size much smaller than the burner diameter. 

On the basis of the supposition that combustion within the primary 
zone was entirely mixing-controlled, the major product species within 
the individual fluid elements were assumed to be present in their 
equilibrium proportions which, for a fixed pressure, inlet temperature, 
and given fuel type, were thus uniquely determined by the local value of 
the fuel: air equivalence ratio and gas temperature. The composition 
of the product gases at any given downstream cross-section of the 
burner could therefore be described in terms of a suitably chosen 
equivalence ratio distirbution function, f(<i>,t), such that f(i}),t)d4) 
represented that fraction of the product gases with fuelrair 
equivalence ratio between 4> and <t> + d<t), with t being equal to the 
transit time of the plug-flow between the primary zone and the 
downstream cross-section. The mean concentration of any major 
component of the product gases, e.g. CO, O 2 , etc., at a given down- 
stream position in the burner was then expressed in the form: 

<[C0]> = /^[C0](T,(j))f{(J>,t)d(() ( 1 ) 

where T is the local temperature— equal to the adiabatic flame 
temperature for no heat transfer. 

It is generally agreed that, for fuels which do not contain 
organically bound nitrogen, the bulk of burner generated nitric oxide 
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is formed in the hot product gases via the well known Zel'dovich 
mechanism. This is an overall endothermic process and, therefore, 
relatively slow by comparison with the fast hydrocarbon fuel 
oxidation chemistry. Provided that the local nitric oxide 
concentration is everywhere less than the corresponding equilibrium 
concentration, it can be shown that the local rate of formation of 
nitric oxide by the Zel'dovich mechanism is, to a close approximation, 
only a function of the local temperature, the local oxygen atom 
concentration, and the molecular nitrogen concentration. By 
assuming that the oxygen atom concentration in the product gases 
was the local equilibrium value and that the local molecular nitrogen 
mass fraction was essentially unchanged by chemical reactions, it 
then followed that the local rate of formation of nitric oxide within 
the fluid elements of the product gases was simply a function of their 
local temperature and equivalence ratio. Thus, the mean kinetically 
limited rate of nitric oxide formation at any downstream cross- 
section of the burner could be expressed in the form: 

<d[N0]/dt> = /^[N0](T,4))f((j),t)d(() (2) 

and the mean cross-sectional average concentration obtained by 
integrating over the transit time: 

<[N0]> = /^[N0](T,4.)f(()>,t)d<|)dt (3) 

In order to apply Eqs. (1) and (3) to calculate the mean 
concentrations,<[CO]> and <[N0]>, a suitable functional form for 
f((j).t) must be chosen. In the modelling studies referred to above 
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the form generally adopted has been a Gaussian distribution 
function: 

f((|>,t) = CVv^ a(t)]exp[-((j) -<<})>) ^/2 a^(t)] (4) 

where <4i> is the mean operating equivalence ratio for the burner 
and 

a^(t) = <(^ - (5) 

is the variance. Pompei and Heywood determined the variation of 
a (t) empirically by assuming the Gaussian form for f((|),t) and 
measuring the mean molecular oxygen concentration along the burner 
length when it was operated at the overall stoichiometric condition. 
At this condition, <[02]> decreased from easily measurable values 
at the upstream poorly mixed end of the burner to zero at the down- 
stream end where turbulence and diffusive mixing had rendered the 
product gases more uniformly mixed. The variation of a (t) was then 
estimated by matching the measured values of <[ 02 ]> with those 
calculated using the expression: 

<[02]> = /^(,[02]{T.<j)) f(<j>.t)d(j) (6) 

By arguing that their mixing parameter, defined to be 

s(t) = o(t)/<4>>, (7) 

was dependent only on the manner in which the fuel and air were 
injected into the burner and independent of <<j)X Pompei and Heywood 
were able to show that this simple plug-flow mixing model was capable 
of describing the effects of nonuniformities in the product gases on 
the burner generated NO and CO concentrations, and that the mixing 
parameter, s(t), could be related, at least empirically, to the 
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overall burner operating variables such as the fuel atomizing air 
pressure. The application of this type of model to investigate 
the effects of changing certain of the operational and design 
parameters controlling nitric oxide emissions from gas turbine 
combustors has recently been reviewed by Heywood and Mikus [5]. 

In this paper we shall address some of the more obvious 
limitations of this simple mixing model, ultimately directing our 
attention to the role of turbulent mixing on burner generated nitric 
oxide formed from organically bound fuel nitrogen for which the 
simple mixing model is inadequate. 

Numerous experimental investigations carried out in both pre- 
mixed flames [6, 7] and oil fired burners [8-12] have shown that 
a large fraction of organically bound fuel nitrogen can readily be 
converted to nitric oxide during combustion. In another paper [13] 
we examined plausible kinetic models which were shown to be quite 
successful in explaining the measured nitric oxide yields obtained in 
premixed flame experiments [6, 7]. On the basis of these models it 
appears that the yield of nitric oxide by combustion could at least 
be minimized if the bulk of the fuel could be burned at fuel -rich 
conditions because, under such conditions, the combustion product 
species which contain single nitrogen atoms have a greater chance 
of being removed from the system by reactions of the type: 

NO + N ->■ N 2 +0 

NH + NO ^ N^O + H 

Indeed, such a view is in keeping with the observations that, whereas, 
those operational methods which seek to minimize burner generated NO 
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by reducing peak combustion temperatures (using such techniques 
as exhaust gas recirculation) are effective for NO formed by 
thermal fixation of molecular nitrogen (the ZeTdovich mechanism), 
only variants on the technique of staged combustion (initial burning 
of the fuel at rich conditions) are effective in reducing the NO 
formed from fuel nitrogen [12], With this in mind, it is clearly 
apparent that in order to construct theoretical models of burner 
flows which incorporate the combined influence of the fuel nitrogen 
conversion kinetics and turbulent mixing, it will be essential that: 

(i) they provide an accurate statistical description of the 
composition nonuniformities in the upstream region of the 
burner where the bulk of the fuel is contained in kinetically 
important fuel -rich fluid elements, and 

(ii) allow for the fact that the net rate of removal of NO 
from the fuel -rich fluid elements is explicitly dependent on 
the concentrations of product species containing single 
nitrogen atoms. 

To briefly amplify on the first point, it is apparent that even 
assuming the strictly one-dimensional plug-flow model, the distri- 
bution of fuel and air cannot be Gaussian initially but must, 
more nearly, be bi -modal with one peak corresponding to pure air 
(cf) = 0) and another corresponding to pure fuel (cj) = °°). The 
question thus arises as to how valid the assumption of a Gaussian 
distribution in equivalence ratios is as used in the modelling 
studies referred to previously. 

With regard to the second point, the reason why Eq. (2) 
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may be used to calculate the kinetically limited mean rate of formation 
of NO is that, according to the Zel'dovich mechanism, the local rate 
of formation, [NO], is explicitly independent of both [NO] and [N] for 
[NO] « [N0]g, i.e., it has a zero order dependence on the concen- 
trations of kinetically rate limited species. When [NO] is not only 
dependent on <{> and T but, also, on the concentrations of rate 
limited species, which is the case for the conversion kinetics 
of fuel bound nitrogen [13], the averaging procedure expressed by 
Eq. (2) is obviously inadequate. 

In an attempt to take a more realistic account of the 
above points, we have employed a stochastic fluid particle 
interaction model of the turbulent mixing process. The model is 
developed from the fluid particle collision model of turbulent 
mixing proposed by Curl [14] to describe mixing in dispersed phase 
systems and, also, by Evangelista et al . [15, 16] and others [17-19] 
to describe mixing in both dispersed and continuum systems. As 
with the simple mixing model, the stochastic model makes no attempt 
to describe the detailed dynamics of the turbulent flow but, 
instead, relies on a composition distribution function to describe 
the nonuniformities in the flow that is itself dependent on an 
empirically determined mixing intensity parameter. However, 
although the method is essentially a numerical one, employing a 
Monte Carlo computational technique, it does allow the shape of the 
composition distribution function to be determined as it evolves / 

with time, given an initial distribution of fuel and air. At the 
same time, it allows one to take account of any general chemically 
rate limited process. Consequently, we have been able to explore 
the validity of the assumed Gaussian distribution function of the 
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simple mixing model and, further, have indicated how a more realistic 
treatment may be made of entrainment or secondary air addition in, 
for example, gas turbine combustors [5, 20]. 

Before developing the stochastic model of turbulent mixing, we 
shall first summarize our experimental measurements of burner 
generated nitric oxide derived from organically bound fuel nitrogen 
[11, 21]. Comparisons between the experimental results and model 
predictions will then follow. 

ir. EXPERIMENTAL MEASUREMENTS OF BURNER GENERATED NO 
FROM FUEL BOUND NITROGEN 

Only a brief description of the burner will be given here since 
it has been described in detail elsewhere [3, 11]. The burner is 
circular in cross-section, 4 in. internal diameter, and 24 in. long. 
Combustion air was supplied at atmospheric pressure and admitted at 
one end through 45° blade-angle swirl vanes. Kerosene doped with 
small percentages of nitrogen containing additives (pyridine, C^H^N, 
and pyrrole, C^H^) , burned at a rate of about 8 lb. hr. ^ was 
supplied through a pneumatic atomizer at one end on the burner 
axis. The air pressure drop across the fuel atomizer was controlled 
independently of the fuel flow rate, and the air flow through the 
atomizer only amounted to about 4 percent of the total combustion 
air. The burner walls were refractory lined to reduce heat losses, 
and the product gases were sampled with the aid of a fully 
traversable stainless-steel probe. The composition of the sampled 
gases was then measured by passing them through an analyzer train. 

All concentration measurements reported here are cross-sectional 


averages. 
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Two series of experiments were carried out in which the total 
air flow was held fixed at 125 lb. hr. In the first series, 
the fuel flow rate was held constant and the fuel atomizing air 
pressure was varied to change the fuel -air mixing process within 
the burner. Low atomizing pressures (p == 12 psig) produced poor 
initial mixing of the fuel and the air and a highly luminous flame 
region was observed; higher atomizing pressures (p - 30 psig) 
produced more intense initial mixing and a blue nonluminous flame 
region. Figure 1 shows the axial profiles of the nitric oxide 
concentration obtained at the two extreme atomizing air pressures 
for the stoichiometric combustion of kerosene doped with 0.51 
percent nitrogen by weight in the form of pyridine. Such results 
are representati ve of those obtained using different fuel flow 
rates and different concentrations of fuel nitrogen, i.e., the 
nitric oxide concentration was observed to rise rapidly to a 
plateau which, for the largest atomizing pressure and the smallest 
nitrogen additive concentration (0.25 per by weight), indicated 
about one hundred percent conversion of the fuel nitrogen to nitric 
oxide; for lower atomizing pressures, the total conversion of fuel 
nitrogen was observed to be considerably less. When approximate 
account was taken of the nitric oxide formed by the Zel’dovich 
mechanism, this was estimated on the basis of previous measurements 
[3] performed at the same operating conditions with undoped fuel, 
the results shown in Figure 1 indicate that about 80 percent of 
the fuel nitrogen was converted to NO at an atomizing pressure of 
29 psig, and at an atomizing pressure of 12 psig the measured NO 
represents about 50 percent of the fuel nitrogen. 
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It is apparent, on the basis of the type of results 
illustrated in Figure 1 that, in addition to the kinetic mechanisms 
responsible for the conversion of fuel nitrogen to NO, the nature 
of the nonuniform turbulent mixing processes within the burner 
must play a very important role in the overall production of NO. 

In the second series of experiments, NO concentrations were 
measured in the exhaust gases, i.e., in the plateau regions of 
Figure 1, over a wide range of fuel flow rates at different atomizing 
pressures using both pyridine and pyrrole as additives. These 
results are plotted in Figures 2a-2c as a function of the operating 
fuel: air equivalence ratio of the burner; each set of data corresponds 
to a different proportion of nitrogen in the fuel. 

First, it is apparent that, as found in previous premixed 
flame studies, e.g., Fenimore [6], the percentage conversion of 
fuel nitrogen to NO is independent of the form in which the nitrogen 
was added to the fuel. The full lines in Figures 2a-2c represent 
the NO concentration which would have been found in the exhaust 
gases for 100 percent conversion of the fuel nitrogen only. The 
dashed lines represent the sum of 100 percent conversion of the 
fuel nitrogen plus that, estimated on the basis of previous measure- 
ments [3], formed by thermal fixation of N 2 . For well-mixed (high 
atomizing pressure) lean combustion (<fi < 1) and the smallest percentage 
of fuel nitrogen additive (0.25 percent nitrogen by weight), the 
measured exhaust NO concentration is seen in Figure 2a to be equal 
to that estimated on the basis of 100 percent conversion of the 
fuel nitrogen plus that obtained by thermal fixation. For larger 
additive concentrations. Figures 2b-2c, the maximum conversion is 
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somewhat less: about 90 percent for the data in Figure 2b, and 

80 percent in Figure 2c. For fuel -rich combustion, ^ > 1, the 
percentage conversion of the fuel nitrogen is seen to fall quite 
drastically. Also, as noted previously, reducing the atomizing 
pressure at a fixed equivalence ratio has the effect of reducing 
the NO concentrations in the exhaust. 

We have previously reported [11] our attempts to measure the 
exhaust concentrations of molecules containing CN radicals; our 
purpose was to establish the possible alternative chemical forms 
of that fraction of the fuel nitrogen which was not converted to 
either NO or N 2 when the burner was operated at low fuel atomizing 
pressures. Briefly, exhaust samples were passed through a gas 
bubbler containing a basic solution of potassium hydroxide in an 
attempt to dissolve such molecules as HCN. Titration of this 
solution with silver nitrate [22] was then followed by the standard 
confirmatory test for ion CN [23]. Analysis of a sample taken when 
the burner was operated at the overall stoichiometric condition, 

<f) = 1 , using fuel doped with 0.51 percent nitrogen as pyridine and 
an atomizing pressure of 12 psig, indicated an exhaust concentration 
equivalent to about 30 ppm of HCN. However, it was shown that 
the efficiency of the system for the dissolution of CN containing species 
low concentrations was poor, so this figure could only be considered 
as a lower limit. We did not pursue this phase of the investigation 
further because, as will be made apparent in later sections, we 
now believe that the major portion of the fuel nitrogen in our 
burner experiments was converted to either NO or N 2 - 
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III. MODELS OF TURBULENT MIXING WITH CHEMICAL REACTIONS 
Corrsin's Model 

The mechanics of the interaction between turbulent mixing and 
chemical reactions has been studied by Corrsin [24-26] for a 
restricted class of flows, i.e., homogeneous, isotropic, constant 
density flows. Although Corrsin's model is not applicable to flows 
of variable density, the results of the simple model will provide 
a useful test of the stochastic models to be presented in the 
following section. 

Following Corrsin, the flow conservation equation for a chemical 
species i with a local concentration (mole cm~^) may be written: 

ar./3t + 9(r.u.)/3x. d r. + R.(r) (8) 

* ■ J o * • 

where x. are the three cartesian position coordinates, 

vJ 

— ~ ^2’ ^i ’ concentration vector of 

the chemical species which compose the flow, is the source term 

representing the production of species i due to chemical reactions, 
and D is the molecular diffusivity which, in Eq. (8) is assumed to 
be the same for all species. 

For a homogeneous flow, i.e., a flow in which the spatial 
gradients of mean properties are identically zero (3< >/9x- e 0), 
the equation for the time rate of change of the mean species con- 
centrations is obtained by taking an ensemble average of Eq. (8) 
over the entire fluid mass: 

d<r.>/dt = <R^(2)> 


(9) 
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The simplest measure of the nonuniformities in composition of such 
a flow is the concentration variance, <y^>» where the equation for 
the time rate of change of the variance is obtained by separating 
T- in Eq. (8) into its mean and fluctuating components, i.e., 

r, " > 

subtracting Eq. (9) from Eq. (8), multiplying through by ^nd 
again taking an ensemble average; this yields: 

d<Y?>/dt = -2D<{3y./3xJ(3Y,-/3x.)> + <Y^R,(I)> 

By further restricting consideration to flows which are isotropic, 
the three mean square derivatives in Eq. (10) may be characterized 
in terms of a single local length scale, s,, analogous to Taylor's 
"microscale" of turbulence, to give 

d<Y->/dt = -(12D//)<Yi> + <YiR^(r)> (11) 

The quantity 12D/)l^, which has the units of inverse time, is a measure 
of the turbulent mixing intensity and will be designated 3(t). 

The simple mixing model, referred to in the introduction, used 
Eq. (9) as the basis for calculating the mean rate of nitric oxide 
formation by assuming that the local composition of the flow could 
be represented in terms of a suitably chosen distribution function 
of equivalence ratios, i.e., f((|),t), and that, according to the 
Zel'dovich mechanism, the chemical source term, Rj^q, was only a 
function of the local equivalence ratio and temperature viz. 

RpjQ = [N0](T,c|>) 


(12) 
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The variation of f(tf),t), or more precisely the equivalence ratio 
variance, along the burner was assumed to be governed by an 
equation similar in form to Eq. (11). However, the assumption that 
3 was independent of the composition was made out of necessity 
and justifieda posteriorie, since it was evident that chemically 
reactive turbulent flows with heat release could not be regarded as 
constant density flows. 

By considering a single species or group of species whose 
concentration, denoted by rf, is unaffected by chemical reactions, 
for example: 

T' = 2[02l + CO] + eCCO^] + [CO] + CH2O] + etc, 

i.e., any chemical element of the reactive mixture, Eq. (11) may be 
written in the simplified form: 

d<yf^>/dt = (^3) 

Therefore, if<y^^> can be measured along the length of the burner, 
the local variation of 6 can be determined from Eq. (13). This is 
essentially the technique used by Pompei and Heywood except that they 
assumed that B could be determined from an equation, analogous to 
Eq. (13), for the variance of the equivalence ratio, viz. 

do^(t)/dt = - 6 a^(t) (H) 

By assuming that the equivalence ratio distribution function was 
Gaussian, see Eq. (4), and that all the major species, specifically 
the oxygen containing species, were present in their equilibrium 
proportions at the local temperature T, they were able to estimate B, 
using Eq. (I 4 ), from measurements of the mean oxygen concentration 
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obtained along the burner length. It is to be noted that the initial 
local temperatures in the combustion zone were taken to be the local 
adiabatic flame temperatures and that these were assumed to decrease 
linearly along the burner so that the total heat transfer rate from 
the burner as a whole matched that estimated on the basis of measured 
outside wall temperatures. Pompei and Heywood thus demonstrated that 
the square of their mixing parameter, s (t) =* a (t)/<(})> , decreased 
in an exponential-1 ike manner, in accordance with Eq. (14)., with a 
high value of 3 in the combustion zone and a lower value thereafter, 
and furthermore, was only dependent on the fuel atomizing pressure. 

A Stochastic Mixing Mode], 

A stochastic model of the mixing process may be developed by 
representing the composition of the fluid at any position in a 
reactor as an ensemble of fluid particles which are identified 
in terms of their individual thermodynamic states. The particles are 
assumed to be sufficiently small with respect to the microscale, l 
(see Eq. 11), that the state of each particle may be assumed to be 
uniform throughout its volume at any instant of time, but the number 
of particles in a given state may vary with time as a result of 
pressure changes, chemical reactions, heat transfer, and mixing 
between particles. 

The thermodynamic state of a representative fluid particle is 
defined in terms of certain of its intensive state properties such 
as the pressure p, the specific enthalpy h, and the species mass 
fractions Y = (Y^ , 'i^, ...» Yj, where m is the total number of species 
contained in the particle. Other representations of the thermodynamic 
state are clearly possible--see, for example, Galant and Appleton [27]. 
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However, the above representation will be used here and in vector form 
this may be written: 

Z = (Y, P, h) 

For a homogeneous , isotropic flow, the state of the ensemble 
can be described in terms of a time dependent distribution function 
of thermodynamic states, f(Z,t). This distribution function is 
simply the density of representative fluid particles in the (m + 2)- 
dimensional thermodynamic state hyperspace of the ensemble, normalized 

in such a way that 

mi-2 

r f(Z,t)dJl = 1 : (dJl = n dZ.) 
n -1=1 

and thus, if we consider equal mass fluid elements,f(Z,t)dfl, is the 
mass fraction of fluid having a thermodynamic state in the interval 
dsi centered about Z. The mean value of any property g(Z), for the 
ensemble is then defined as: 

<g> = /n g(I)f(I.t)da 

In the absence of mixing, temporal changes in f(^,t) may be formal- 
ly described by the continuity equation for the concentration of fluid 
particles in the thermodynamic state space: 

m+2 

3f(Z,t)/3t + E aCf(Z,t)Rf(Z)]/3Z. = 0 (15) 

i=l 

The second term may be expanded to explicitly describe the changes in 
f(Z,t) caused by pressure variations, heat transfer (other than con- 
duction between fluid elements, which will be treated with species mix- 
ing), and chemical reactions, respectively, viz. 
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m+2 

Z 3[f(Z,t) R.{Z)]/8Z. = 3[f(Y,p,h;t) dh/dt]/9h 

i=l ^ 

+ 3 [f(l.P.h;t) dp/dt]/3p 

m 

+ Z 3[f(Y,p,h;t) R.(Z)]/3Y. 
i=l ^ ^ 

For an adiabatic {no heat transfer to the reactor walls), constant 
pressure flow of the type which we shall consider, the first two terms 
of the above expansion are obviously equal to zero. 

The change in the distribution function due to turbulent mixing may 
be considered by incorporating a mixing term into Eq. (15) to create a 
master equation of the form:. 

3f(Z,t)/3t + 9[f(Z,t) R^(Z)]/9Z. = Of(Z,t)/3t}^.^ (16) 

Any mixing term which conserves the total number of points in the thermo- 
dynamic state space, i.e., conserves mass, may, in principle, be applied 
in Eq. (16). In the following development the mixing term will be 
evaluated using a model which describes the mixing process as sequence 
of interactions between randomly chosen groups of n representative 
fluid particles which mix completely with one another to form n new 
particles having identical thermodynamic states. Between interactions 
chemical reactions are assumed to proceed within each fluid particle 
as though it were a closed system. Thus, a class of master equations 
is defined in which the mixing term is derived in much the same way as 
the binary collision term (n = 2) of the kinetic theory Boltzmann 
equation in which the particles are assumed to be statistically inde- 
pendent of one another between collisions. 

If the ensemble is assumed to be composed of equal mass fluid par- 
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tides, a mixing interaction between n particles with thermodynamic 
states described by the vectors will yield n 

particles with the same thermodynamic state, viz. 

Z, = (z/') + z/^) + ... . Z>>)/n, i . 1,2 m + 2 (17) 

where the pressure, p = is assumed to be constant throughout the 

flow field. Then the probability that any random interaction between 
n particles in states ..., Z^”^ will produce n new particles 

in the particular state ^ is: 

n f(Z^^^t) f(Z^^\t) ... f(Z^^\t) 

X n 6{{Z.^^^ + z.^^^ + ... + Z>^)/n - Z.} da^^^ ... 

1=1 1 T 1 1 

where the Dirac delta function, is equal to unity when Eq. (17) 

is satisfied and zero otherwise. Similarly the probability that mixing 
will remove a fluid particle in state ^is: 

n f(Z,t) f(Z^^^t) ... f(Z^*^"^^t) 

X (1 - 6{Z. - + ... + Z.^""^^)/(n-l)>] ... 

i=l ' ' 

If the frequency at which a fluid particle is mixed completely with 
(n-1) fluid particles, which may be a function of the local thermodynamic 
state of the gas, is denoted by co, the net rate of change of f(^,t) due 
to mixing is obtained by multiplying the probabilities that mixing will 
increase or decrease f(^,t) by w and then summing over all possible 
interactions. Thus, Eq. (16) becomes: 
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m+2 

3f(Z,t)/3t + E 3[f(Z,t) R.. (Z)]/3Z. = 

7 = 1 

^ ... f(Z^"U) 


X ["^n^ 6{(z/^^ + ... + z/"h/n- z.}]d$^^^^ ... dQ 
1=1 ^ ^ 


(n) 




fi 


( 1 ) 


r(n-l) 


■ ^ (1) •'* ■'^.(n-i)2“ f(r .t) ... f(r \t) 


X [\ (1 - 6{Z. - + ... + z/""^h/n})] ... (18) 

i=l ^ ^ ^ 

For the moment we shall limit our consideration to constant density 
flows in order to show that the results derived from Eq. (18) are 
formally quite consistent with Corrsin's, i.e., Eqs. (9) and (11). 
Assuming constant density, p, the simple transformation 

r. = pY./m^ , 

where m^ is the molecular weight of species i, converts the mass fraction 

-3 

of Eq. (18) to the concentration, e.g., mole cm , used by Corrsin. 

Noting the assumption that the molecular diffusi vities of all species are 
the same and that a Lewis number of unity is implicit to our interaction 
model, the mixing frequency, w, may be assumed to be independent of 
the fluid particle state. An equation for the time rate of change of 
the mean thermodynamic state is then obtained by multiplying Eq. (18) 
by ^ and integrating over all This yields 

d <Z.>/dt -<R.(Z)>= 0 (19) 

which is Corrsin's result. 

The equation for the variance of state property Z^. is obtained by 
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multiplying Eq. (18) by = (Z. - <Z->)^ and again integrating over 
all Z, viz. 

d<z.^>/dt - <z^.R^(Z)> = -a)(n-lKz^^> (20) 

By comparing Eq. (20) and Eq. (11) we see that they are the same if the 
mixing frequency, w(t), is related to 6(t) through the expression: 

U)(t) = B(t)/(n-l) (21) 

Thus, in the constant density limit the fluid particle interaction model 
yields the same result as Corrsin's more conventional approach. 

If the initial distribution function is known, Eq. (18) may be 
solved numerically using well known Monte Carlo methods [28] for any 
general chemical reaction kinetics. We shall assume that, in applying the 
stochastic mixing model to our burner, the burner may be modelled as 
a simple plug flow reactor as before. The mixing frequency o) is assumed 
to be equal to 6/(n-l) and to be independent of the local gas composition. 
Thus, we regard w as being empirically determined in the same way as it 
was when the simple mixing model was applied. The numerical Monte Carlo 
integration procedure for Eq. (18) is illustrated in Fig. 3 and is briefly 
described in the following paragraphs. 

First, an ensemble of N equal mass fluid particles is chosen to repre- 
sent the incoming fluid composition of the burner. Each particle is as- 
sumed to be of uniform composition. A fraction of these particles, (N^/N), 
is taken to represent the flow through the fuel atomizer and to have a 
composition determined by the equivalence ratio of the atomizer jet 
(<|) »1 ) . The remaining particles, (N - Ng) , represent the pure combustion 
air, = 0. Each particle is assigned a number index, 1 s i s N, and a 
pseudorandom number generator is then used to generate n different indices. 
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i„, .... i {i. f i.), thereby identifying the particles to be 
I c n j K 

involved in the first interaction. They are allowed to mix completely, 
separate, and react according to the appropriate chemical rate equations. 
The mixing interactions are computed sequentially, therefore, this 
first interaction represents a forward step in time equal to 

At = l/w(t)N (22) 

where (jj(t) is the empirically determined, time dependent mixing frequency. 
The chemical rate equations are then integrated over the time interval 
At. A new group of n fluid particles is then chosen at random from the 
entire ensemble of N fluid particles to be involved in the next mixing 
interaction; in this way the calculation proceeds through the burner. 

At any time during the mixing process, the mean composition and 
other mean properties of the fluid in the burner may be evaluated by 
taking an ensemble average over the N fluid particles, e.g. 

Z = (1/N) 2 (23) 

i=l 

Furthermore, the shape of the evolving distribution function, f (^,t) , 
may be determined by plotting the ensemble population as a function of 
any one of the thermodynamic state variables or any other local 
property . 

Equation (21), i.e., 
w(t) = g(t)/(n-l) 

guarantees exact correspondence between Corrsin.'s description of the 
turbulent mixing process and that based on the stochastic n-body 
interaction model for the case of constant density flows. We shall thus 
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proceed on the assumption that in applying the stochastic model 
of mixing to the burner flows, the calculations can be performed 
using the simplest model, i.e., n = 2, the binary mixing model, 
and delay numerical justification until the end of the following 
section where comparisons between the results will be made using 
n = 2, ... .6. 
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IV. STOCHASTIC MIXING MODEL RESULTS 
The Evolution of the Distribution Function 

For a burner which admits the fuel and the air separately and 
which is to be treated as a truly one dimensional plug-flow, the 
initial distribution function for the fuel and the air will, ideally, 
be bimodal rather than Gaussian as assumed in the simple mixing model. 
Nevertheless, the assumption of a Gaussian distribution function did 
result in good estimates of both the NO and CO emissions [3]. We 
shall now examine the reasons for the apparent success of the Gaussian 
distribution function assumption using results obtained from the 
stochastic mixing model. 

The calculations were performed for the case of stoichiometric 
combustion using the computational procedure outlined in Figure 3, An 
ensemble of 1000 fluid particles was used to represent the flow. 

Figure 4 illustrates the shape of the distribution function computed 
at four downstream positions within the burner as a function of the 
dimensionless time: 

e = /6(T)dT = C(t)/N (24) 

0 

where c(t)is the total number of two-particle mixing interactions 
which have taken place within the ensemble prior to time t. At early 
times, when the fluid particles have only mixed an average of once, 
e = 1, most of the fluid is still contained in either very fuel rich 
or very lean fluid particles with only a small fraction at or near 
the stoichiometric condition. Only after each fluid element has mixed 
three or four times does the distribution function approach a Gaussian 
distribution, as indicated by the dashed lines in Figure 4. 

Now according to the Zel'dovich mechanism, the maximum rate of 
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formation of nitric oxide occurs only within a narrow range of equivalence 
ratios about the stoichiometric condition, i.e., 0.8 5 tf) ^ 1.2, where 
the adiabatic flame temperature is highest. The mass fraction of the 
fluid which, at any time, is responsible for the major portion of the 
NO production is, therefore, 

1.2 

/ f(4»,0)d<|) 

(|i-0.8 

In Figure 5 we have compared this fraction, as calculated using the 
simple mixing model (full line) and the stochastic model (point values), 
as a function of the dimensionless time 9. The agreement between the 
two is obviously very good even at early time and thus, it is apparent 
that the reason for the success of the simple mixing model is not that 
the Gaussian distribution function has the correct overall shape, but 
that it correctly describes the fraction of the fluid which is contained 
within the important narrow range of equivalence ratios about the 
stoichiometric condition. 

The simple mixing model also provided reasonable estimates of 
the mean carbon monoxide concentration within the burner. By assuming 
that the local CO concentration was the local thermodynamic equilibrium 
value, it is obvious that only the fraction of fluid which was fuel 
rich could contain significant CO concentrations. Thus, in Figure 5 
we have also compared the mass fraction of the fluid with ^ > 1.05, 

using both the simple mixing model (broken line) and the stochastic 
model (point values). For this case it is seen that only for times 

9 > 3, do the two predictions agree; for 6 < 3, the simple model 
overestimates the fraction of fuel-rich fluid and may thus overestimate 
the mean CO concentration in the poorly mixed upstream region of the 
burner. By contrast to the case of NO, the model prediction of the 
mean CO concentration is dependent on the shape of the distribution 
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function. To determine the actual extent of this region of discrepancy 
between the two model predictions requires that we know 6 as a function 
of position within the burner. 

Now complete combustion can only occur after the fuel and the air 
have mixed to some degree, and the range of cj) within which complete 
combustion is likely to occur is not well defined. However, based 
upon estimated flammability limits for kerosene [29], we have assumed 
that complete combustion occurs if 0.3 < 4> f 3.0, and that no combustion 
occurs outside these limits. Accordingly, Figure 6 shows the mean 
oxygen concentration, calculated using the stochastic model for the 
case of stoichiometric combustion, as a function of the dimensionless 
mixing time, 9. By comparing the values of the mean oxygen concen- 
trations measured as a function of position within the burner [3] with 
those contained in Figure 6, we were thus able to empirically estimate 
the variation of 6 within the burner for different fuel atomizing 
pressures as illustrated in Figure 7. It can be seen that the variation 
of 9 implies that B(t), see Eq.(24), decreases from a higher value 
in the upstream region of the burner to a lower one downstream, in 
accordance with Pompei and Heywood's observations and, furthermore, it 
is apparent that even at the lowest fuel atomizing pressure, 9 approaches 
a value of 3 within the first diameter of the burner length. Thus, 
the predictions of the mean CO concentrations obtained from the simple 
mixing model should be reasonable over the major downstream portion of 
the burner. 

In order to incorporate the actual kinetic model for NO production 
(the Zel'dovich mechanism) into the stochastic model of the burner 
flow, it would have been necessary to have taken account of the heat 
losses from the burner. Indeed, the levelling-off of the mean NO 
concentration with axial distance within the burner, which was observed 
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experimentally [3], was directly attributable to the fact that the 
highly temperature sensitive rate of NO formation decreased with 
axial distance due to the accumulated effect of heat transfer. The 
approximate iterative method used by Pompei and Heywood to account 
for heat losses is not so easily incorporated into our stochastic 
model. For the burner under consideration, the primary mode of 
heat transfer from the combustion gases to the burner walls was 
conduction. Although, in principle, this could be modelled by a 
consideration of the fluid particle interactions with the wall, 
it would have involved lengthy iterative procedure to match the 
calculated heat transfer rate to the total heat transfer rate which 
Pompei and Heywood estimated by a consideration of the natural con- 
vection and radiation at the outside wall of the burner, see Section IV 
reference [3]. Furthermore, without more detailed measurements of the 
heat transfer rate as a function of axial position, it would be necessary 
to introduce rather arbitrary and, possibly, questionable assumptions 
regarding the degree of temperature accommodation which could be 
achieved between the fluid particles and the wall during interactions. 

For these reasons and because we have shown that the simple model 
provides a satisfactory description of the composition nonuniformities 
over the major portion of the burner length, we have not attempted to 
incorporate the Zel'dovich mechanism of NO formation into the stochastic 
model of the burner flow. A more satisfactory test of this aspect of 
the problem is provided by modelling the NO formation in a more 
adiabatic burner, i.e., a burner without significant heat transfer 
effects, such as the gas turbine combustor recently considered by Mikus 
and Heywood [20]. We shall briefly discuss their results in the following 
section. 
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Gas Turbine Combustor 

Nitric oxide formation in conventional gas turbine combustors may 
be modelled ty assuming the primary zone to be a partially stirred 
reactor to which the simple mixing model, as described earlier, may 
be applied [1-3,5]. Very little nitric oxide is formed outside the primary 
zone due to the rapid quenching which occurs with the admission of the 
secondary air. Mikus and Heywood [5] showed that extensions of this 
basic model could not adequately describe the NO emissions over the 
full range of operating conditions of a new combustor design concept 
presently being developed by NASA [30]. In this design, the NASA 
modular swirl can combustor, fuel and part of the combustion air (about 
10 percent) is introduced through 120 small swirl cans arranged in three 
concentric rings at the same axial station; most of the remaining 
combustion air is introduced as flow around the swirl cans. 

At very fuel-lean operating conditions, combustion occurs primarily 
in the small recirculation zones in the wakes of the individual cans. 

Each reaction zone can thus be treated separately as a partially stirred 
reactor using the simple mixing model, but including the reverse of 
the NO formation reactions. By suitably choosing the mixing parameter, 
s (see Eq. 7), the NO emissions can be correlated as shown by the 
lower broken line in Figure 8 [5]. As the overall fuel:air equivalence 
ratio approaches unity, the individual reaction zones merge so that com- 
bustion may then be viewed as occurring throughout the entire volume of 
the combustion chamber. The correlation illustrated by the upper broken 
line in Figure 8 was obtained by treating the entire combustor as a 
partially stirred reactor. It is evident that extrapolations based on 
these two limiting versions of the simple model do not correlate the NO 
emissions at intermediate operating conditions. 
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Mikus and Heywood [20] have recently applied our stochastic 
mixing model to describe the flow in the NASA swirl can combustor 
in an attempt to correlate the NO emissions over the full range of 
operating conditions. The composition of the gases leaving the 
recirculation zones of the individual swirl cans was assumed to be 
described by the simple mixing model as applied to a partially stirred 
reactor. Mixing between this gas and the secondary air which flowed 
around the swirl cans was then described using the stochastic mixing 
model. The actual rate of entrainment of the air by the product gases 
was not treated separately. Instead, the function 3(t) was chosen so 
that the calculated NO emissions gave a best fit to the measured values. 
Used in this way, 6(t) is a combined measure of both the turbulent 
mixing intensity and the entrainment rate. The full line in Figure 8 
illustrates the extremely good correlation which was achieved. A more 
detailed discussion of these results is contained in the recent paper 
by Mikus and Heywood [20]. 

Burner Generated Nitric Oxide from Fuel Nitrogen 

In a previous paper [13] we proposed a simple kinetic model to 
describe the yields of nitric oxide measured in the product gases of 
premixed fuel-rich flames to which a variety of organic fuel nitrogen 
compounds had been added [6,7]. According to the model the fuel 
nitrogen is first assumed to be rapidly distributed in its local 
equilibrium proportions at the adiabatic flame temperature among the 
species which contain a single nitrogen atom, i.e., NO, N, NH, NH 2 , etc.- 
these we designate collectively as RN. The most abundant RN species 
for local fueliair equivalence ratios less than about 1.8, is NO. If 
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the amount of organic nitrogen in the fuel yields a "constrained" 
equilibrium concentration [27] of RN in excess of the thermodynamic 
equilibrium value, then its removal will proceed by reactions which 


ultimately produce N 2 , e.g. 




N + NO 

N 2 

+ 0 

(1) 

NH + NO 

N 2 O 

+ H 

(2) 

N + NH -> 

N 2 

+ H 

(3) 


and, accordingly, the rate of removal will be described by a rate 
equation of the form: 

d[RN]/dt = -2k(cl>,T}[RN]^ (25) 

where the effective rate coefficient is 

k(<t.,T) = k^([N]/[RN])g([NO]/[RN])g + .... 

+ k3([N]/[RN])g([NH]/[RN])g (26) 

and the suffix e denotes equilibrium concentration ratios. 

For the simple plug-flow burner under consideration in which the 
fuel and the air are admitted separately, it is to be anticipated 
that most of the fuel will first be burnt in fuel-rich fluid elements. 
Thus, in applying the stochastic model of the mixing process to the 
burner, we see that the fractional reduction in the concentration of 
RN in a fuel -rich fluid particle is given by an expression obtained by 
integrating Eq. 25 with respect to time, viz. 

[RN]/[RN]^^^q = l/d+At/x) (27) 
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where At Is the time measured between interactions and the characteristic 
reaction time, t, is defined to be; 

T = V(2k((l,,T)[RN]^^^Q) (28) 

Figure 9 illustrates the variation of the effective rate coefficient, 
2k((}),T}, as a function of the fueltair equivalence ratio for the 
adiabatic combustion of kerosene fuel at atmospheric pressure. The 
rate coefficients: k^, k^, etc., which are contained within the 

definition of k{<|),T) were evaluated using the expressions given in 
reference [13]. 

Figure 10 illustrates the variation the NO concentration, calculated 
on the basis of c ompl ete thermochemi cal equilibrium, as a function of 
fuel:air equivalence ratio — also for the case of adiabatic kerosene 
fuel combustion. The broken line illustrates the nitric oxide 
concentration in the product gases which could be formed if the fuel 
nitrogen (0.51 percent by weight as N in the fuel) were converted to 
the species RN on the basis of the constrained equilibrium assumption. 
According to the kinetic model, the rich fluid particles which fall 
within the cross-hatched area of Figure 10 will lose fuel nitrogen at 
a rate given by Eq. 25. However, after mixing with air or more fuel- 
lean fluid, the composition will change and the representative particles 
may then be located at positions to the left of the cross-hatched 
region. In this region the NO concentration will be less than the 
corresponding equilibrium concentration (indicated by the full line) and 
thus the rate of removal of the RN species will be exceeded by the 
formation of NO via the Zel 'dovich mechanism. As we have indicated 
previously, significant rates of NO formation by thermal fixation only 
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occur for fluid with fuel: air equivalence ratios between about 0.8 and 
1.2. At lower equivalence ratios the rate of formation is very slow 
due to the low local temperatures. For that portion of the fluid which, 
by further mixing, is characterized by fuelrair equivalence ratios less 
than about 0.5, the local NO concentrations may again be far in excess 
of the local equilibrium values. However, due to the exceedingly low 
concentrations of all of the RN species other than NO at such low 
equivalence ratios, further removal of NO via reactions (1) - (3) will 
be virtually halted on time scales comparable with the flow transit 
times through the burner. 

With this qualitative description of the mixing process in mind we 
see that, to a first approximation, the rate of formation (destruction) 
of NO within the burner may be modelled by first assuming that the NO 
derived from the fuel nitrogen is formed at a rate independent of the 
amount formed by thermal fixation. Thus, the fuel derived NO may be cal- 
culated using the stochastic model of the mixing process by applying 
Eq. 26 to describe the rate of RN removal within each fluid particle 
between interactions. In carrying out this calculation, no account need 
be taken of the effects of accumulated heat transfer from the burner 
which gives rise to an exhaust gas temperature drop of about 200°-300° K 
[3]. This assumption is justified by observing that the rate of NO 
removal as described by Eq. 25 is relatively insensitive to such small 
temperature variations — certainly much less so than is the rate of NO 
formation by the Zel'dovich mechanism. The total nitric oxide concen- 
trations within the burner may then be obtained by adding to these results 
the additional NO which is formed by thermal fixation— the latter being 
estimated on the basis of Pompei and Heywood's simple mixing model 
calculations which include the effects of heat transfer. 
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Obviously, this combined model for calculating the total NO 
formation will provide an upper bound to the NO concentrations within 
the burner. The reason for this is that the model takes no account of 
the fact that some fraction of the thermally derived NO, which is formed 
in fluid with a near stoichiometric composition, will be mixed into 
fuel -rich fluid and may, therefore, be removed by reactions (1) - (3) 
at a rate described by Eq. 25. The magnitude of this error may be 
examined by considering a fuel-rich fluid particle which has a residence 
time At between mixing interactions and an initial single nitrogen species 
concentration: [RN] = [RN]^ + [RN]^, where [RN]^ refers to the con- 

centration derived from the fuel nitrogen and [RN]^ refers to the con- 
centration formed by prior thermal fixation. According to the combined 
model, the RN concentration after time At will be 

[RN]^^. = [RN]^ + [RNyd + 2Atk{(j>,T)[RN]f) (29) 

whereas, in reality, the removal reactions should also affect the NO 
formed by thermal fixation, viz. 

[RNl^^ = ([RN3t+ [RN]f)/(l + 2Atk((t>,T){[RN]^ + [RN]^)) (30) 

In the limit, At/x = 2Atk(<|),T)[RN] « 1, the two results are obviously 

identical. However, for At/x » 1 , we see that 

[RN]'j = [RN]j*[RN]^^ (31) 

i.e., the difference between [RN]^^ and [RN]^^ may amount to as much 
as [RN]^. Thus, we are in a position to calculate both an upp*er bound 
to the NO emissions, i.e., the result of the combined model, and a lower 
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bound — obtained by neglecting the contribution due to thermal fixation. 

Calculations were carried out with the stochastic mixing model 
and the fuel nitrogen conversion kinetics model, as described above, for 
the burner operating conditions that corresponded to the experiments 
reported in Section II. Values of the mixing intensity, B, were 
deduced by matching the calculated mean axial oxygen concentration 
profiles with those measured by Pompei and Heywood. Satisfactory fits 
to the data at each atomizing pressure were obtained by choosing an 
initial value of the mixing intensity, S-j , which remained constant in 
the upstream portion of the burner for a mean transit time t^. There- 
after, the mixing intensity assumed a lower value & 2 > which appeared to 
be independent of the atomizing pressure. These empirically determined 
values of B-| » given in Table I as a function of the 

fuel inlet air atomizing pressure. 

Axial profiles of the calculated mean NO concentrations for the 
case of overall stoichiometric combustion (0.51 percent nitrogen by 
weight in the fuel) at the two extremes of the fuel atomizing pressure 
are compared with the experimental measurements in Figure 11. The 

full lines, estimated on the basis of the stochastic model, illustrate 
the contribution of the fuel nitrogen alone; the broken lines add to 
this the NO formed by thermal fixation and were estimated using the 
simple mixing model. For the well mixed case (p = 29 psig) the broken 
line closely approximates to the measured profile. The reason for this 
agreement is that at high fuel atomizing pressures the residence time 
of the fuel in the rich regions of the flow is curtailed by rapid 
mixing, thereby, minimizing the effect of reactions (1) to (3) which 
reduce the fuel nitrogen to N 2 * By contrast, it can be seen that at 
the lowest atomizing pressure consideration of only the fuel nitrogen 
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kinetics, as shown by the full line in Figure 11, yields good agree- 
ment with the measured NO concentration in the exhaust. The reason 
for this stems from the fact that mixing takes place sufficiently 
slowly at the lowest atomizing pressure that a large proportion of 
the RN species (about 40% of the original fuel nitrogen) can be 
converted to in the fuel-rich fluid. Clearly, this case approxi- 
mates the asymptotic limit for which the residence times of 
representative fuel -rich fluid particles are much longer than the 
characteristic chemical conversion times. 

One feature of the calculated axial NO profile which is not 
reflected in the measurements is the maximum which occurs at about 
one diameter downstream from the fuel inlet nozzle. This discrepancy 
suggests that the one dimensional plug-flow model does not accurately 
represent the position-time history of the flow in the upstream region 
of the burner. We can only speculate that recirculation, which is 
necessary to stabilize combustion in the actual burner, is responsible. 
Certainly, a recirculation zone would have the effect of smoothing 
out the peak in the measured cross-sectional average NO concentration. 

Figures 12(a) - 12(d) show comparisons between the calculated 
and measured mean NO concentrations in the burner exhaust for the full 
range of fuel atomizing pressures as a function of the overall fuel: air 
equivalence ratio. For the case of well mixed combustion (p = 29 psi ) 
the prediction of the combined model (broken line in Figure 12(a)) is 
clearly in good agreement with the measurements over a large proportion 
of the range of equivalence ratios investigated. Again, this indicates 
that the intense mixing at this condition does not allow sufficient 
time for the fuel nitrogen to be reduced to N 2 within the fuel-rich fluid 
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elements. For the lowest atomizing pressure (p = 12 psig) it appears 
that the contribution of thermal fixation to the NO exhaust emissions 
is small at all equivalence ratios since the measurements fall 

closely about the stochastic mixing model prediction which only accounts 
for the fuel nitrogen conversion kinetics. By including the comparisons 
which have been made in Figures 12(b) - 12(c) for the two intermediate 
fuel atomizing pressures, we recognize an obvious trend in the functional 
dependence of the exhaust NO concentrations on both the overall 
equivalence ratio and the fuel atomizing pressure, namely: the equi- 

valence ratio at which the combined model ceases to correlate the 
measured NO emissions decreases as the fuel atomizing is lowered. 

Even at the highest atomizing pressure, Figure 12(a), the two data 
points at equivalence ratios greater than 1,3 are seen to lie closer 
to the lower bound estimate. In order to calculate the exhaust NO 
emissions which fall midway between the upper and lower bound estimates 
it would be necessary to include thermal fixation in the stochastic model 
and, although this is easily incorporated, it would require that an 
accurate account of the heat losses from the burner be included. 

Figures 13(a) - 13(b) and 14(a) - 14(b) show comparisons between 
our upper and lower bound model predictions and the measurements of the 
exhaust NO concentrations as a function of equivalence ratio for fuel 
nitrogen additive concentrations of 0.25 and 1.06 percent by weight 
nitrogen, respectively. The same general trends are observed for the 
well mixed and poorly mixed conditions as were described above. However, 
it is now apparent that the lower bound prediction overestimates the NO 
yield by about 20 percent for the case of the highest fuel nitrogen 
additive concentration and the lowest atomizing pressure. We can only 
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speculate that the reason for this stems from the inadequacy of our 
simplified model of the fuel nitrogen conversion kinetics. The 
assumption that the fuel nitrogen is initially distributed amongst all 
the possible single nitrogen species, including N and NO, is probably 
suspect forthe very fuel-rich fluid elements— a point which we discussed 
in our previous paper [13]. Until a better and more exact understanding 
of the fuel nitrogen conversion kinetics has been developed, we cannot 

profitably pursue this application of the burner mixing model further. 

Approximate numerical justification for the validity of the 

relationship 

o)(t) = 3(t)/(n-l) 

for variable density flows with finite rate chemistry and heat release 
is provided by the results shown in Figure 15a and 15b. For each 

mixing model, i.e., n = 2, 3, 4, 5 and 6, values of , t^ and 
were chosen so that the computed values of the mean molecular oxygen 
concentration were closely matched to the measured values as a function 
of axial position within the burner. By comparing these values, see 
legend to Figure 15, it can be seen that in no case do the values 
differ by more than 12 percent; this difference is probably within 
the statistical uncertainty of the calculations caused by the finite 
number of particles in the ensemble (N = 500 for all results shown 
in Figure 15) . 

The close agreement between the mean NO concentration profiles 
shown in Figure 15b implies that each model is capable of correctly 
representing the mixing-compos’it’on-t'sme history of the flow as a 
whole. In view of this it seems that the stochastic nature of the 
turbulent mixing process is the central most important aspect of 
the mixing models and that the details of the interaction, i.e., the 
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number of fluid particles involved or the equal mass assumption, 
is relatively unimportant. Indeed, it is likely that other stochastic 
model representations of the mixing term in the master equation, 

Eq. 16, could be made to work equally well. However, from a practical 
point of view it appears that the simple binary mixing model is 
justified in describing the type of burner flow considered here. 

V. CONCLUSIONS 

In this investigation we have developed a stochastic model of 
turbulent mixing for a plug flow reactor in which mixing is represented 
by random fluid particle interactions. The model was developed for 
the purpose of investigating the role of composition nonuniformities on 
burner emission characteristics. As with previous mixing models such 
as the partially stirred or plug flow reactor concepts, no attempt is 
made to describe the detailed dynamics of the turbulent flow. Instead, 
the composition of the flow as assumed to be described in terms of a 
distribution function of thermodynamic states which evolves with time 
by the action of general rate limited chemical processes and mixing. The 
mixing process is characterized in terms of an empirically determined 
"mixing frequency" which has been shown to bear a direct relationship to 
the more familiar mixing rate intensity as used, for example, by Corrsin. 
The model thus represents the flow as an ensemble of perfectly stirred 
reactors, each characterized by its own thermodynami c state, which are 
created and destroyed by fluid particle mixing. Between interactions, 
chemical reactions are assumed to proceed at the finite rates dictated by 
appropriate kinetic models as though each reactor (fluid particle) were 
a closed system. 
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On the basis of this stochastic mixing model we have shown 
that the simple plug flow mixing model (developed by previous authors 
to investigate the role of turbulent mixing on burner generated thermal 
NO and CO), which assumes that the composition nonuniformities can be 
described by a Gaussian distribution function in the local fuel: air 
equivalence ratio, is justified. 

The stochastic model has also been extended to include the 
combined effects of turbulent mixing and secondary air entrainment on 
thermally generated NO in gas turbine combustors. The model enabled 
measured NO emissions of a new modular combustor design to be 
correlated over a wide range of fuel: air ratios in terms of a single 
empirically determined mixing/entrainment rate parameter. 

Finally, we have used the stochastic mixing model to predict 
rate limited upper and lower bound estimates on the yield of burner 
generated NO formed both by thermal fixation and from organically bound 
fuel nitrogen. The resulting comparison which we were able to make 
between the model predictions and the measured NO yields served to 
illustrate that the stochastic mixing model can be used to relate burner 
emissions to such operating variables as fuel atomizing pressure, overall 
equivalence ratio, and fuel nitrogen additive concentration. Further- 
more, we conclude that quantitative agreement of the model predictions 
and the measurements in this latter application was only limited by 
insufficient information on the local rate of heat loss from the burner 
and by possible inadequacies in the kinetic model of fuel nitrogen 
oxidation and reduction. 



39 - 


ACKNOWLEDGMENT 


This work was supported in part by the National Aeronautics 
and Space Administration under Grant NCR 22-009-378 and by the 
Environmental Protection Agency under Grant R-800-729-03-01 . Part 
of the work was presented at the A.I.Ch.E. 66th Annual Meeting under 
the title: Nitric Oxide Formation from Fuel Bound Nitrogen in an 

Oil Fired Burner. 

We would like to express our thanks to several of our colleagues 
for helpful discussions, notably Professor D. P. Moult, who suggested 
extending our investigation of the stochastic mixing model beyond the 
two- particle limit. 

One of us (R.C.F.) wishes to acknowledge partial financial 
support from M.I.T. for a Sloan Research Traineeship. 



- 


TABLE I 

Mixing Intensity 


P (psig) 

3^{sec"^) 

tg(sec) 

B2(sec 

29 

400 

.0205 

60 

20 

300 

.021 

60 

17.5 

252 

.0216 

60 

15 

200 

.023 

60 

12 

112 

.025 

60 
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FI6URE CAPTIONS 


Fig. 1 
Fig. 2 

Fig. 3 

Fig. 4 
Fig. 5 

Fig. 6 
Fig. 7 
Fig. 8 

Fig. 9 

Fig. 10 
Fig. n 


Axial variation of NO concentration (ppm, wet basis) 
during combustion of kerosene with 0.51 percent nitrogen 
by weight added as pyridine. 

Exhaust NO concentrations for combustion of nitrogen 
containing fuels. Filled symbols: pyridine additive; 

Open symbols: pyrrole additive. 0> P ■ 29 psig; 

O, p = 20 psig;^, p = 17.5 psig;A, p= 12 psig. 


Computer code flow diagram for stochastic mixing model. 


Note: 


(i) 


where ZjJ 


(i) 


IS 


intensive state property of i^ fluid particle. 


Distribution function f(<l)) as a function of dimensionless 
time 9. Broken line - Gaussian distribution function. 


Mass fraction of fluid within specified ranges of 
equivalence ratio plotted as a function of the dimension- 
less time 0. Points represent the stochastic mixing 
model evaluation and the lines represent the Gaussian 
distribution function evaluation. 

Variation of the oxygen concentration with the dimensionless 
time 6 calculated for stoichiometric combustion using the 
stochastic mixing model. 

Variation of 9 with axial distance along burner evaluated 
using oxygen concentration measurements of Pompei and 
Heywood [3]. 

Comparison of measured and calculated nitric oxide emissions 
from the NASA swirl can combustor [20] at an operating 
pressure of 5 to 6 atm. Broken lines: simple mixing 

model calculations;. full line: stochastic mixing model 

calculation. 

Variation of effective rate coefficient for RN removal 
reactions as a function of fuel; air equivalence ratio for 
adiabatic combustion of kerosene. 


Equilibrium nitric oxide concentration as a function of 
equivalence ratio for adiabatic combustion of kerosene. 

Comparison of measured and calculated axial nitric oxide 
concentration profiles for combustion of kerosene with 
0.51 percent nitrogen added as pyridine. Full lines: 
fuel nitrogen contribution alone; broken lines: combined 

model predictions. 





Fig. 12 Comparison of measured and calculated NO emissions 

from combustion of kerosene with 0.51 percent nitrogen 
added as pyridine, #, and pyrrole, O. Full lines; 
fuel nitrogen contribution by stochastic model; broken 
lines; combined model predictions. 

Fig. 13 Comparison of measured and calculated NO emissions from 
combustion of kerosene with 0.25 percent nitrogen added 
as pyridine, #, and pyrrole, O* Full lines: fuel 

nitrogen contribution by stochastic model; broken lines: 
combined model prediction. 

Fig. 14 Comparison of measured and calculated NO emissions from 
combustion of kerosene with 1.06 percent nitrogen added 
as pyridine,# , and pyrrole, O. Full lines; fuel 
nitrogen contribution by stochastic model; broken lines: 
combined model prediction. 

Fig. 15 Comparison of stochastic model calculations for n = 2, 3,... 6. 

(A) 0 Measured oxygen concentrations. 

Calculated oxygen concentrations: 
q : n = 2, = 200 sec-1, t = 0.024 sec, Bg = 60 sec'l ; 

□ : n = 3, B-| = 200 sec~\ t = .021 sec, ~ 

n = 4, = 200 sec'\ x = .021 sec, 63 = 60 sec“^; 

A‘* n = 5, a^ = 200 sec"\ x = .021 sec, ^2 ~ 60 sec'^ ; 

V: n = 6, a^ = 175 sec~\ x = .023 sec, 62 = 62.5 sec'\ 

(B) Calculated fuel nitrogen concentrations: [RN], open 

symbols; [NO], filled symbols. 
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Figure 3* Computer Code Flow Diagram for Stochastic Mixing Model 


\ . . . ) , where ^ 
is the intensive state property of the i^^ fluid 


Note: Z{^) = 


particle. 
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